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ABSTRACT 

Many algorithms for polynomial least squares approximation of real-valued function on a 
real interval determine polynomials that are orthogonal with respect to a suitable inner prod- 
uct defined on this interval. Analogously, it is convenient to compute Szego polynomials, i.e., 
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he appropriate to deterinine Szego polynomials in algorithms for least squares approximation 
of real-valued periodic functions by trigonometric polynomials. This paper is concerned with 
Szego polynomials that are defined by a discrete inner product on the unit circle. We present 
a scheme for downdating the Szego polynomials and given least squares approximant when a 
node is deleted from the inner product. Our scheme uses the QR algorithm for unitary upper 
Ilessenberg matrices. We describe a data-fitting application that illustrates how our scheme can 
be combined with the fast Fourier transform algorithm when the given nodes are not equidistant. 
Application to sliding windows is discussed also. 
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Downdating of Szego Polynomials 
and Data Fitting Applications 

CI.S. Ainnioi’ WM3. Gra<^i>^ L. Reichel" 



Abstract 

Many algoritlims t'or polynomial least squares approximation of a real-valued function 
on a real interval dcteririinc i)olynomials that are orthogonal with respect to a suitable inner 
product defined on this interval. Analogously, it is convenient to compute Szego polynomials, 
i.e,. polynomials that are orthogonal with respect to an inner product on the unit circle, 
when appro.ximating a complex-valued function on the unit circle in the least squares sense. 
It may aUo be appropriate to determine Szego polynomials in algorithms for least squares 
approximation of real-\ allied periodic functions by trigonometric polynomials. This paper 
is concerned wnh hzego polynomials that are defined by a discrete inner product on the 
unit circle. \\« jireseiit a scheme lor downdating the Szego [lolynomials and given h.'ast 
Mjuares appro.Miiiaiit when a node is dcdeied from the inner product. Our scheme uses the 
C)ll algoriihm lor uuilar\ upper lle^enbeii; matrices. We describe a data-fitting application 
that illustratc'N how our scheme can be combiinM with ih* *' fast Fourier transform algorithm 
when the giuMi nodes arc not equidistant. Application to sliding windows is discussed also. 



1 Introduction 

Lot !)(' a of distinct nodes on thr unit circle and h't be a set of positive 

weights. Introduce for complex- valued functions fj and h defined .ii the nodes the discrete inner 
[iroduct on the unit circle 

<n 

((j.h'lni ■■= ( 1 . 1 ) 

k=[ 

whore the bar denotes complex conjugation, rolvnomials that are orthogonal with respect to 
an inner product on the unit ciiclc are known as Surjo polynomials. Let denote the 

family of ortbonnrmal Szej>6 polyiioiuials with respect to the inner product (1.1), where 6j is of 
degree j and ha:> a positive leading coelficiont. Fhe polynomials <i>j satisfy the Szego recurrence 
rrlations 



Oi)[z) = Oo{z)=\/rr^, 

* l)c*partmeul of .M ;ii hemalical Sciences. ,\’oi tliciii lHiiiois I niversity, DcKalb, IL 60115. 

’ Depai 1 ineiil of .\Li i lieiii.il ics. .Vaval Fo<i gindiia i o School. Monicicy, ( \\ 9 3910. 

* Department ol Mai heiualics and Coiupnier Science. Kent Male University, Kent, OH •112<12. Research 
>iipported bv a .Naiional Research Council fellowsliii^ and NSF grant DMS-9002884. 
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rr_,+,Oj+i(-) = ro_,(r) + -j + iOj(-), 0<;<m-l. (1.2) 

^T^+iOj + i(c) = r-:,+,<pj(r) Oj(;). 

wlioi'c the recurrence coefficieuts 'ij+\ € C :ii\d + i > 0 are deteriniued by 

/A •V'"' 

^0 = ^\j = I 2^ (f’L j 

7 ,. I - ( 1 . 3 ) 

^^ + 1 = [l-hv+il'j • 0 < j < m. 

See, lor example, Gioiiaiuler and Sze.12,6 [ 4 . (.'liaj)tcr 2]. It can be shown by induction that the 
auxiliary polynomials Oj in (1.2) satisiy 

Oj(z) = z^Oji 1/c), 0 < j < m , 

whoie 6 j{z) denoto.s the [)olynoinial obtained by conjugating the coefficients of 6;(r) in the 
power basis. Since the measure on the unit circle that defines (1.1) has precisely rn points of 
increase, we have |7^| < 1 for i < j < m and |-„i| = 1. The coefficients 77 are Iciiown as Scliur 
parameters, and we refer to the as the associated complementary parameters . Although the 
com[)lcmentary parameters are mathematically redundant, we retain them during; computations 
ill order to avoid luiinerical instability. Note from ( 1 . 3 ) that (75 is the total weight of the measure 
that defines (1.1), and that is the leading coelficient of d)j for 0 < j < m. 

Tor later reference we aUo define the polynomial 







mC^m— 1 ( ■“ ) 


(1.4) 


riieii 










= fi(- 

A = l 


: - -A.), 


(1.5) 


In particular, (<Am. d>j) 


= 0 for 0 < j < m. 






Example 1.1. Let zi; 


:= exp(2a/(A: - l)/m) and u’ 


i ■= 1 •‘or 


1 < k < m. where i := \/— T. 


riieu ao = aj = 


1 and 7j = 0 for i < j < m. 


nud frn = 


1. Thus. d>j{z) = for 


1) < j < 7n. and dmC-) 


= 77 |-^/-(r^^ - 1). 




□ 



ll.'-yllm := 

•) 



Introduce the discrete norm 



and let doiiolo the sol of nil polynoini.ii.^ oi (l('«j;rcc nt most n - 1. Let / be a »iven 

complex- valued function defined at the nodes r; . and consider the problem of computing the 
polvnomial p^i-i € Hn-i* fot* >ome n < ni, that ■^aiishes 

\\f - l>n-l\\m = i>V‘ ll/-/>llm- (1-6) 

The solution Pn-i of the minimization problem ( l .(>) can be expressed conveniently in terms of 
the Szego polynomials Oj as follows. Introduce ilie \ ecior 

f = . (1.7) 

and the w x // malnx () =■ 

(ji,j o,_i i ^ < J < (i.S) 

u here := n r’. .Vote that the matrix () Imn oi ( lioiiorinal cobunns. i.e.. 0‘Q = /. where 
(y := () ^ . Lot the vc'Ctor a = ^o i . o _> n, l>o dolim'd bv 

a := (rf . ( 1 . 9 ) 

riien the i)ol\'iiomial ]>n-i be written a.s 

V, . (1.10) 

= I 

where the rourier coellicients (\j are indepeiuk'iii of n. 

It is the purpose of tliis i)aper to present an <dgorithm for downdatuu] the reciiirence coef- 
licienls 7^ <ind a,. <\^ well the Fouiier c(Kdll^i('nt.^ n . when an node- weight pair is removed 
from the inner product (1.1), Our algorithm is ba>ed on 1 he observation that t he coluinns of the 
mntiix 0 are eigenvectors of unitarv llesscuilHuu niairix defined by the recurience relations 
{ L.‘ 2 )-( 1.3). This makes it i)ossible to downdaie iIk' coeincients 77. Cj and Oj by applying the 
i}[l algorithm for unitary llessenbeig maliice>. poNcmU'd in [ 5 ], with the node to be removed 
as shift. Details are described in Section 2. 

We remark tliat the problem of uiidahiuj iln‘ ( oellicTmts 7^, Gj and aj when a node- weight 
pair {-m-hi ^ added to the inner [)ioduci ( 1 . 1 ) is discussed in [S]. The updating problem 

can l)e sol\*ed by using an inverse QR algorithm lor unitary llessenberg matrices; see [S, 1]. 

.\ssume that we wish to determine the polvnomial /Ui_i. given by (l.O)-(l.lO) with n < m 
when the set of m nodes in the inner piodiict (I.l) is a subset of the set of .V equidistant 
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nodes {exp( 2 t/j/.V with k := .V - ni > 0 small. The weights are all assumed to be 

unity. Then it may l)e attractive to compute by hrsi coinj)uting the polynomial interpolant 
/>.v-j € n.v_i oil the sot of A’ (Miiruli.^uiiit llO(l('^ hy ii>iiig the fast Fourier transform (FFT) 
algorithm, and then deteriuiiiiug Pm-i from /^v_i \>v <ipj)lying our downdatiiig scheme. The 
Fourier coefRcieius of the polyaoinial aj)[)ro.\imaiU /^i_i are then equal to the first u Fourier 
coefTicients of Pm-i- This approach can similarly be applied to trigonometric polynomials. 
Details are described in Section 3. Computed e.Kamples are presented in Section 4 and Section 5 
contains a summary. 

Updating and downdating of polynomials ap[)roximants pn-i when all the nodes are 
real has received considerable attention in the literal ure; see Scott and Scott [9] and references 
there. A collection of algorithms for updating and downdating based on orthogonal polynomials 
is presented by Flhay et al. [3]. Algorithms for updating are also considered in [G, 7]. 

2 Downdating of Szegd polynomials 

The connection between the Szegd polynomiaU determined by (1.1) and a unitary llesseii- 
berg matrix can be seen as follows. Using the basis of Szegd polynomials, we can write 

I < J < m. (2.1) 

;=i 

where + lor 1 < j < ///. and + = I - b<'t ///, ; '.= 0 for 3 < j-h2 < k < m, and define 

the upper Hessenberg matrix II = define i he unitary matrix U - i by 

//A-., )/r/. . (2.2) 

Substitution of c = zi;, I < k < m, into (2.1) yiedds the e(piation 

IlT'u'^ = r'r\. (2.3) 

or, equivalently, 

II = U'\U. (2.4) 

where A := diag[ri.c_> Thus. ![ is a uiiiiary u[)per Hessenberg matrix with positive 

subdiagonal elements that is uiiicpielx* delenniin'd by the inner ])roduct (1.1). Also note that 
the first n columns of U make up tlu' matrix () chdined bv ( 1.8). .Mgorithms for the solution of 
the least-squares problem ( l.(i) ran t li(M('fore \)o vi('W('(| in ((nnis of the spectral decomposition 
( 2 . 1 ). 
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It is fairly siraiiilit forward to show, by iisiiiy the recurrence relations ( 1 .2 )-( 1.3), that 11 
can be written as a product of //? olcMiicntary unitary transformations that are defined by the 
recursion coelRci(Uit.s , and a, for the o. . Wo have 

7/ = t / I ( * 1 j 0 2 ( V2 ) * * ‘ f - 1 i 7m - I ) m ( 7m ) ? (2.5) 

where the 6 'j( 7 ;). 1 < 7 < <tre m x m Givens matrices 

C-. 

and G'|^( 7 „J i.^ tiu' diagonal inairix 

^ ■■= diay[i . i 

We refer to the rc'pre.^oiitation (2.5) a.^ the Schiir parametric form of //. 

The development of efficient algorithms for oigonproblems for unitary Hessenberg matrices is 
facilitated by the fact that e\ory unitarv nj)per Hessenberg matrix with nonnegative subdiagonal 
elements has a unlipie Sdiiir pa rametei ization. l or example, when the implic.tly shifted QR 
algorithm is aj)plied to find the spectral decoinposii ion of a unitary Hessenberg matrix, a 
sec[uence of intermediate nnitaiy l[e.''S(Miberg inatiice> 1 .^ generated that converges to a diagonal 
matrix. In [o] the QR algorithm for unitary 1 [e.s.-onbcrg matrices is formulated in terms of the 
Schur parameters of the iiitcrmedi.ue iiiatrice^. This results in an implementation that requires 
only 0(m) arithmetic operations ])cr iteration on a matrix of order m. 

Assume for the moment that the matrix II and scaling factor ctq are given, but that the 
nodes Zk nnd weights icj of the inner product ( 1 . 1 ) are not explicitly known. Then it follows 
from (2.4) and (2.2) that the nodes arc the ('igcnvalucs of //, while the normalized weight 
Wk!<yo is equal to the first component of the normalized eigenvector corresponding with Zk* (We 
a.^sulne that each <‘iuciivoctor is scaled \n liav<' unit Kuclidoan length niul a noiuicgaiivc first 
com])onciU.) The unitary llosscnberg algorithm can be used to compute the matrix A and 
vector i'e\, without compuiiiig the of l\ wiih iioi more than 0(in) arithmetic oi)crations 
p(M‘ iteration, 'riiroiigliout this papoi cj (h'iioi('> i lu' axl^ x'cetor in C'^7 I'lie QR algorithm 
determines one node- weight pair at a time, ami for each pair computed the order of the unitary 







1 



0 



ii[)pcr llcsseubcrg matrix is reduced by one. so that the reduced Hessenberg matrix corresponds 
with the node-weight pairs that have not yet been determined. 

W^e remark that in the case that tlie nodes of the discrete inner product arc real, tlicn the 
analogue of the matrix // is a real symmetric tridiagonal matrix T with positive subdiagonal 
elements. This matrix T contains recurrence coellicicnrs for orthonormal polynomials that 
satisfy a three-term recurrence relation. Golub and W'elsch [4] proposed the use of the QR 
algorithm for symmetric tridiagonal matrices for the computation of the node-weight pairs 
associated with T. This algorithm also determines one node-weight pair at a time, and reduces 
the order of T when such a pair has been found. 

Conversely, ilie construction of // from the inner product (1.1) can l)c regarded as an 

inverse eigenvalue problem. In particular, given the matrix A = diag[ci.z> and the 

vector qo := we can perform a seciuence of elementary unitary similarity 

transformations whose composition results in an m x /n unitary matrix U such that the matrix 
on the right-hand of 



■ 1 0'^‘ ■ 








■ 1 o'’’ 




ic 




.0 C*. 




.Mu A . 




.0 / 






II 



is of upper Hessenberg form with positi\e .Mibdiagonal elcmieiits. (The ★ represents an arbitrary 
scalar that remains iiiiclianged. ) In other words. C‘qo = ^r^tl U'AU is r. unitary upper 

Hessenberg matrix II, Conseciiiently. II is the inaiiix corresponding with the imter product 
(IT). 

The transformation of A to II can be p<M formed using 0{m?) arithmetic operations with 
the inverse unitary Qll (lUQR) algoiiihin described in [l]. The lUQR algorithm is an updating 
procedure because it incorporates node-weight pairs one at a time. After the stage of the 
algorithm has been executed, the / x / unitary Hessenberg matrix corresponding with the inner 
})rodiict determined by ilie first j nodes and weights has been obtained. 

In [8] the approximation problem (l.O) is solved using the lUQR algorithm to construct 
the Schur parameters of the unitary Hos.senberg matrix II . The elementary unitary matrices 
are accumulated against the vector f during the algorithm to obtain the vector Fourier 
coefficients a = //’f without explicitly Ibriniiig I’, In this way we obtain the interpolating 
polynomial Pm-i lhat is the solution of (1.6) with u = m in 0(m") arithmetic operations. Of 
course, the partial sums of the Fourier expansion of the interpolating polynomial yields the 
.solution (1.10) of (l.G) for each n < in. Moreover, if one is interested in computing only Pn-i 



lor n < then the al^oritlim can 1)C curUnlcd mj that only Oirnn) aritluuetic operations are 
required for tlie coin])iitalion of the parametoi- and See [<'^] for 

dor ails. 

Xow assume that we have >oKed rlie least-^(|u«u ('> piol)leiu ( 1.6) with // = m hy the luctliod 
described in [sl. >o that sots ol' Schiir paraiiuM fu >. and of corn j)leinent<uy parameters 

{^;}yLo tiorrespoiuliiie; with the inner [Product ( 1 . 1 ). as u ell as sets of Fourier coolficieuts 
of the interpolatiim' polynomial Pm-i are explicitly known. In the downdatin^; problem, we seek 
to ^()lve the corrcs[)oii(liue; least-Mpiare^ problem with one term deleted from (1.1). In particular, 
let f be the node that is to bo deleted from the inner product (1.1) and let w bo the s(|iiare*root 
of the nssociat('d W('i.L;ht. In order to 'Impliry ^oim* lormiilas that follow, we assume without 
loss of geiicralitv that i = and ic = a*,,,. 

Introduce the inner product 

- 1 

1 ^/.// ) , _i := ^ ' :. l/i )//’^ (‘2. 7) 

and discrete norm 

\\‘o seek the polynomial /ha-> € 11,/i-j "uch that 

1 / ~ /^n - j i — 1 — min _/ — /yjrn-i • (•^•b) 

/'*5l 

Denote the fainilv ol orthonormal S/<'eo pol\ iiomlaU with positive leading cocincicmt associated 
with (2.7) hy .\1no lot and (hmotc the sets of recurrence coeiricieiits 

for the Oj. and hu .\ := diag^Cj ' <‘U(1 (},, := , wlune rrj, ;= {a' — 

is the total weight of the moa.Miie that delines t2.7). Then from the above discussion, 
there is a unitpie unitary llcssmiborg matrix II and a unicpic unitary matrix U such that 

iqo = 



and 

II = r'.UX (2.9) 

where denotes the 7''* axis vector in C'^“^ Moreover, the om iirrencc coefficients 'jj and 
cTj arc the Scliiir paramel(‘i'S and comj)lemcnt ary [)aiametcrs. rc>|n ( 1 ivcly, of 7 /. The optimal 
polynomial is them ‘aivem by 

r) = \ niOk-\[Z), 



' = 1 



where the vector of rourier cocKicienis a ,oi.r») is given by a := (/"f and 

f := [wif{zi) U%n-[f(-rn-\ )]^ 

Our scheme for (lowiulatiiig is based on the ob>ei\’ation that II and U can be computed 
from II and U by applying one step of the ()R ali;oiirluu with ‘‘ultimate” shift i. together with 
some permutation similarity transformations, lo determine a unitary matrix H' such that 



T 



1 0 
0 II’- 






(Toej 

II 











L O'^' 
0 li^ 


II 









aoei 

w 



aoef 

II 



w 

0 



( 2 . 10 ) 



Then 



and 



r 0 



a 

wf{z) 



= (;ll- 



= 11 "a . 



( 2 . 11 ) 



Tlius. tlie dowiulated Fourier coefRcients d;. are obtained by accumulating each elementary 
unitary transformation against a. An elRcieiiL implementation is obtained by using the QR 
algorithm for unitary Hessenberg inaliiccs described in [o]. 

We now assume tliat i c/ for some /,!</< nr, and let w := in/. One step of the QR 
algorithm with shift : aj)phod lo the niaiiix II determines a unitary upper Hessenberg matrix 
W such that 



//' := r*//r = 



II 0 



because i is an eigenvalue of II . It is easy to s('c. however, tliat the QR algorithm applied to II 
will not yield the re([uired II , because the vector will not he transformed as required by 
( 2 . 10 ). On the other hand, an HQ algorithm lor //, in which the transforming matrix V would 
be a product of Givens matrices in the reverse order of (2.12) below, would yield the desired 
similarity transformation. 

Instead of modifying the unitary Hessenberg (^U algorithm of [ 5 ] to perform an RQ iteration, 
we apply the QR algorithm to the unitary l(<'ssenberg matrix II ^ := JII^J, where J = 
‘ • *Gi] is the reversal niaf ri.r o\' oviU'V in. (l is (Msily seen that if H is given by ( 2 . 5 ), 

then 

( 7l • 2i 7: ) • • • ^ •' ,n - I ( */ I )^m(7m)i 



where 7^ lm-,ri<n for 1 < 7 < m and •= 1 w. • I'ho application of the QU algorithm on 
//^ is equivalent with that of the RQ algorithm on II . One iteration of the QR algorithm with 



^hift i applied to u^ciicratos a iiiiitarv upper 1 Icsscnbcrg matrix 

1 ^ M ( <'' 1 : i ) ■ ■ ■ ^ - 1 ( (\n - 1 m ( ) 

Mich that 



r = 



//' 0 

0^ . 



;2.i2) 



(2.13) 



Moreover, S,n can be taken to be an arbitrary iiniinodiilar number because deflation has taken 
place [5]. Let := (1 — denote the compleinentary parameter to dj for I < J < rn. 

Observe that only the last two conii)oncnts of Woe^ are nonzero, and that \Sm\ = 1 can 
be chosen so that these two coiuponents are given by 



a — 1 

^lU-l 



rrj - t 





0 




— 1^0 








_ l<^m-lko _ 



Filially, we transforiii i he rigliL-hand >ide ol (2.1.d) by similarity using 



r 



./ 

0^ 



0 

1 



where J is 



the ie\*ersal lualrix of order /n - 1, and irans|>os(' the result. In this wav we transform // to 



'^IIW = 



II 0 



0^ 



. where U' = ./\* 



. Moreo\'Pr, IF’^oei = 



<^ 0 ^ 771 - 1^1 

l^m- 1 1<^0 



, and 



./ 0 
0^' 1 

by the unicpieness of the reduction, II” = //, rror^_i = do, and |f^_i|ao = w. The vector of 
Fourier coefficients a is then determined from (2.11) by applying W* to a incrementally. 

Observe that onr downdating i)rocedure requires knowledge of the node to be deleted but not 
of the corresponding weight. We can therefore compare the computed value of xh to the actual 
\alne <r/ in order to assess the accuracy of the computations. If iv is not close to u’/, then the 
(low nd at ed polynomials and prn-: might not be accurate, .\nother accuracy check is provided 
by the computed value of /Un in the algorithm. This (juantity is the diagonal element of 
the upper triangular matri.x in the ()U factori/ai ion of //^' — i/, which is mathematically zero 
when i is an eigeii\'aliie of //. If the coin[)uied value of Pm is not *‘tiny'\ then the computed 
downdated polynomial might be inaccurate. We iherefore consider Pm ^nd lu as being part of 
the out[)Ut of the algorithm. 



0 



Algorithm 2.1 (Downclating by removing llie node-weight pair { = ,w} ) 
Input: m. -• 

Output: = l - Pint 

r., .im-1 



7,];lv 7,n[7,n-,]r=V; K]7=v := K-.];=V; 

[^*jlj = l [‘^m-j + 1 ]j = l • 



, o 

II 


= I; 






j : = 


1.2... . 


, m — 


I do 


Pj ••= 






+ 7/j 


if;> 


; 2 then | 


:= Tj. 










‘S : = 


(iA,-. 


+ 7j<^j 


-i )!Pj 


Mj : = 




+ r/r 


1 + 1 •' 


‘■‘i + i 




: -f ^ 


/+! ■' 


^ •• = 


(Vi ^ 




-I )IPj 


Ij '■ = 


- 






;= |i^ 


m-l + 


7m ^m — 


il.- 






^m— 1 • — ^m — l/^m> 



h =■ |(^ 



m- 1 






7m_i := 7m-i/|7m-i|.' := O; ( parainctci' coiTcctioii) 






— 7m-l [7m-j- 1 lj = l 






rn — 2 . 

1 • 



= KV„ 



1 m— 1 . 
-j\j=\ ’ 



The algoritlim o\‘or\vrite.s the Fourier coelficioiUs with iutermedinte (pmiititics. Al- 

gorithm 2.1 requires 0{ in) arithmetic operations ( =^. /) and the evaluation of w — 1 square 
roots. 

We have already noted that if the nodes are real, then the analogue of the matrix II 
is a real symmetric tridiagonal matri.x T with positive subdiagonal elements. Similarly with 
Algorithm 2.1, downdatiiig of the ortlioiiormal polynomials associated with T, as well as of 
Fourier coefficients, can be carried out by algorithms based on the QR algorithm for symmetric 
iridiagoiial inairices. This observation may bo new. 

In certain data-fitting applications it may be desirable to update the polynomial Pm-n 
given by (1.9)-(1.10) with n = m, by loplaciiig certain node-weight pairs. This can be carried 
out by successively removing an node-weight priir by .Algorithm 2.1, and then adding a new 
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node- \v('i‘^lit pair ii>in^ uiic .sto[) ol Ali;oiiiliiu d.i ai A]. This coinbinauoii oi al^oriiliins \*lelds 
a slidinij; window schoiiie. 

3 Downdating and the FFT algorithm 

Wdien the nodes are equidistant on the unit circle and all weights U‘1 arc unity, then 
interpolation and least-squares aj)proxiination ol a given function by algebraic or trigonometric 
polynomials can l)e carried out rapidly l>y the FFT algorithm. This section describes in some 
detail how our downdating scheme mav be combined with the FFT algorithm to achieve rapid 
schemes for interpolation when the nodes are (^'^sentially equidistant. More precisely, we will 
consider the case when the set of nodes in the inner product (1.1) is a subset of the 

set of ecpiidistant nodes {exp( ♦ '^'‘dicro k := .V - m is a small positive integer. In 
our Operation count w(' will a-^uine that u is iiidcp<Mi(l('nt of in. The weights trr in (1.1) are all 
assumed to be uiiitw 

Lot / (leiiole .1 t cUii pl< x- X .1 1 iM.'d lillicl'oii. '(>-«• \,iiii<‘N /(C;.), 1 1_ explicitly 

known. We remark iliai a rej)i oseiit a : Ion ol ii,(' int (u pohii ing polynomial /Un-i € IFn-i 
Newton form can be computed in Ol/n*) ariihiiieric opeiations. The Xaiidermoiule solver by 
Bjorck and Pereyra [2] can be used to determine a representation of />,n-i nt terms of the 
monomial basis, and reciuires also 0(m") arithmetic operations. Our scheme only requires 
0(m log 77?) arithmetic operations and yields a representation of Pm-i basis of Szego 

pol\'nomials, 

I.et denote the ((implement of iIk' in ) }^'^_“ and let 

j'(z[ ) := 0, 1 < h: < n. riius. / is deliiied at i h(' A 1 0 ( )i s of u nity exp(2:r/f j- ! )/.V). i < j < .V, 
and we can compute the T'otirier coefhcienis of the polynomial px-\ € that inler[)olates 

/ at these nodes by the FFT .dgorithiu in (A -X log A’) = 0{m log m) aiithmetic operations. 
U.sing the Schui parameters given in Example 1.1. we apply Algorithm 2.1 ac times to eliminate 
the nodes !</«-*< from the inner product. This requires 0(r77) arithmetic operations 
and \uelds the Fourier coeflicients of the desii('(| interpolating polynomial The Fourier 

coeincients of the least squares approximants /a,._i G Un-\ '‘'‘'ilh n < rn are. of course, a subset 
of the Fourier coellicieni s of / a„^_|. 

.A scheme closely relatc'd to the one outlined abo\e can be used to compute trigonometric 
[)olynomials rapidly, l.et r;. and z[ be the imih's introduced above, and dcTine := arg(r/^), 
!</.’< m. and := arg(c[.). 1 < /r < /». .XFo .i.>.^nme that m is odd. i.e.. m =. 2r 4- 1. Let 
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f(0) be a real- valued rniiclion deruied at ilie llod^'^ and let ) •= 0. 1 < h < n. We wish 

to compute a trigonometric [)ol\'uoiuial iiO) G >i)nn{l.cos 0 cos r^.sin 0 sin 7^} that 

minimizes the sum 



kzzl 

This can be accomplished as follows. First solve the minimization problem 



mi 






:V 



-U0-U/Ad2 



J=l 



by the FFT algorithm in 0(w log m) arithmetic operations, and denote the optimal coefficients 
by Cij. Remove the nodes 0[, 1 < A; < k, by ap[)lying the downdating scheme k times to the 
polynomial p 2 r+K(-) •= XTyLu"" where .V - I = ‘2r + k. This rc(|uires 0{m) arithmetic 

operations and yields the polynomial /;o,. ^ flo,.. The desired trigonometric polynomial is then 
given by i(6) := - - <'xp(/^). 



4 Computed examples 

W^e now present the results of some nmuerical experiments with our downdating procedure. 
These experiments were perfoinicd in FORTRA.V on a SparcStation SLC at Northern Illinois 
University, on which there arc a[)proximately 7 and IG significant decimal digits in single- 
precision and doiible-[)iecision calculations, le.spectivcly. 

The first set of ex[)eriinents compares the accuracy of the downdating procedure with that of 
the updating procedure lUQR described in [S]. We ini)iit A' uiiimodular nodes positive 

weights and comi)lex function values For any positive integer Jii < A^ let 

am denote the vector of Fourier coefricients of the solution /^m-i of (i.G) with n = nu and lot 
gm denote the vector of Schur [)arameters determined by the inner product (1.1) and recurrence 
relations (1.2)-( 1.3). 

We first compute \'ectors a.v <'^nd g/v using an im[)lementation of the lUQR algorithm in 
single-precision arithmeiic [s]. We thou repeatedly ap|)ly a single-precision implementation of 
Algorithm 2.1 to coin])iito and g^ for deert'asing values of m. The application of Al- 
gorithm 2.1 removes ilu' node fiom rin' inner |)ifninrt to comi)Ui.o the solution of (I.G) 

with n := m := A^ - /r. After (xich downdatiiig step, wo calculate the relative error in i*e., 
||am - fimILVIlamllz, inul the error in g^, i.e., ||g,H - gMi||- 2 , where ||x ||2 denotes the Euclidean 
norm of x G C'^. We also solve each [)rol)lem of (u'der m using the lUQR algorithm in single- 
precision arithmetic and compute the resulting error's. The results of the lUQR algorithm in 
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tloiible-[)i ccision ariilunelic arc n-ed ns exact answers in the error calculations. The following 
tables displax’ the iosiilling errors for the dowiidatiiig procedure (DD) and the updating proce- 
dure (UD). In all of the cx[)criiaeuts. each function value /(-:;) has its real part and imaginary 
pan randomly generated according to a uniform distribution in [-5,5]. 

\Vc first choose the nodes to be the N := 300 roots of unity Zj := exp(2Tz(j - l)/N), 
I < j < A\ and uniform weights u’j := 1. Table 1 shows the errors for problems of order 
ni := y -k for various values of /c. Table 1 also shows the results with the same choice of nodes 
and weights, except that the nodes are permuted in a random way. This permutation changes 
the nodes that are deleted as well as the order in which the nodes are added in the updating 
procc'dure. It should bo noted I hat the eirors in the downdating procedure can be expected to 
iurrease as k increases. Table 2 shows that ^iluihu• results are obtained with the same set of 
nodes and randoiulv generated weights in (0,10). 

Tiible 3 shows the results obtained with uniform weights and .V := 300 nodes Zj \ — 
<'xp(.T/(j - 1)/.T). 1 < 7 < Ad in [0,~), both in their original order and in a random order. 
Here again, the dowudating procedure perfoiin> well. 

Table 4 shows the results when the initial 300 nodes are randomly selected points on the 
unit circle. In this example, the error in the downdating procedure displays a sudden increase 
at k = 10. In this step the error in the coinputcMl dowiulated weight, \ tc - ^i;|, was greater than 
10~H Ob.serve that the error incurred at l\\\^ dowudating stop propagated to the subsequent 
dowudating steps, but the errors in a.^ and ^('em to grow graduallx' as k increases, i.e., as 
III := X - k decrea>e>. Other experiments with random nodes produced similar results. It 
should be noted that in otir experiments, a large error in the computed weight did not always 
coincide with a large jump in the errors in VL^n and g,„. It is clear that more work needs to be 
done in order to iinder.stand the numerical aspects of the updating and downdating ])roblems. 

Our final experiment tests the accuracy and speed of the procedure described in Section 
3 for dowudating the FFT. The .\'-point FFT is used to obtain the Fourier coefficients of the 
interpolating polynomial pw_i. .A raiulondy selected set of nodes is then removed from the 
inner product using .Alvoriihm 2.1. d'his «'xi>eriiueut was run with a radix-2 FFT subroutine in 
single precision arithiiH'tic with := 102 1 and .V := 20 IS. Table 5 shows the computed error 
after k dowudating steps for various values of k. As above, we use the results of the lUQR 
algorithm in double precision arithmetic as exact answers for error checking. \Vc also display 
the amount of CPU time required by the FFT with k downdating steps and the time required 
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Table 1: 300 nodes with equispaced arguments in [0,2;r); uniform weights. 





nodes in 


order ol increasing; ari;ninent. 




nodes in random order. 






relativo on or m a,., 


error 


1“ grri 1 


1 relative error in 


error 


i“ gm 


k 


OD 


FD 


DD 


ID ii DD 


UD 


DD 


UD 


0 


0.1221::-03 


O.I22E-03 


0.120E-03 


0.12GE-0;! 


0.2S1E-04 


0.'2S4E-01 


0.310E-04 


0.310E-0*! 


10 


0.1 iii:-03 


0.121E-03 


0.295E-03 


0 1 lOE-03 


0.335 E-0-1 


0.290E-01 


0.5 15E-0 4 


0.310E-94 


•20 


O.MlE-03 


0.115E-03 


0.-105E-03 


0 118E-03 


0.173E-04 


0.314E-04 


0.2 11 E-03 


0.3C1E-04 


30 


O.lloE-03 


O.lOOE-03 


0.322E-03 


0. 1.5*1 E-03 


0.550E-04 


0.272E-04 


0.272 E-03 


0.331E-01 


*10 


0.112E-03 


O.OGTE-Ol 


0.310E-03 


0.1.55E-03 


0.943E-04 


0.257E-04 


0.359E-03 


0.374E-04 


50 


0.112E-03 


0.9 10 E- 01 


0 329E-03 


0.164E-03 


0.964E-04 


0.259E-04 


0.3-25E-03 


0.437E-04 


70 


0.124E-03 


0.102E-03 


0.578E-03 


0.171E-03 


0.174E-03 


0.-258E-04 


0.424E-03 


0.422E-04 


90 


0.13-2E-03 


0.105E-03 


0.725E-03 


0.173E-03 


0.251 E-03 


0.235E-04 


0.677E-03 


0.422E-04 


110 


0.151E-03 


0.107E-03 


0.3C5E-03 


0.1 6*1 E-03 


0.375E-03 


0.238E-04 


0.753E-03 


0.383E-04 


1 130 


0.1.51E-03 


0.1 17E-03 


0.-2S9E-03 


0.160E-03 


0.465E-03 


0.226E-04 


0.109E-02 


0.359E-04 


' 150 


O.lGGE-03 


0.127E-03 


0 321E-03 


0.1. 50 E-03 


O.GOOE-03 


0.18GE-01 


0.119E-02 


0.463E-04 



by the single-precision ILT^IL algorithm on the problem of order m = X — L\ It is interesting to 
note that the downdating procedure produces >iil)stanli;dly more accurate answers faster than 
the lUQR algorithiii ev(Mi for moderately .^ized walues of Ic. 

5 Conclusion 

New algorithms are presented for discrete least squares approximation by algebraic poly- 
nomials on the unit circle and by trigonometric i)olynomials. The algorithms downdate the 
approximant as well as recurrence coefficients for Szego polynomials when a node-abscissa pair 
is removed from the inner product. The schemes are based on the QR algorithm for unitary 
Ilessenberg matrices. A particularly attractive application is to combine our algorithm for 
trigonometric polynomials with the fast Fourier transform algorithm. This yields a fast scheme 
for computing trigonometric least squares approxiinants wdien the nodes of the inner product 
form a subset of cHpiidistant nodes on [0,2a). 



14 



Tabic 2; 300 nodes witli c(iuispaccd ar^unioius m [0.2~); random weights in (0, 10) 





nodes in 


order ol iiiero 


asiiii;' ari;umonL. 


nodes in random order. 




relative error in a,n 


error 


HI 


relative error in k\m 


error 


ni gm 


k 


DD 


UD 


DD 


I'D 1 


DD 


VD 


DD 


UD 


u 


0.207E-03 


0.20 7 E-03 


0.21-1E-03 


0 211E-03 . 


0.511E-0-1 


0.511E-04 


0.572E-04 


0.572E-04 


10 


0.206E-03 


O.1S3E-03 


0.33.') E-03 


0,1s iE-o;t 


0..763E-0-1 


0.537E-0-1 


0.814E-04 


0.744E-04 


20 


0.199E-03 


0 171E-0.3 


0.-117E-03 


0.1 92 E-03 


0.G29E-01 


0.512E-01 


0.1 lOE-03 


0.722E-01 


:\0 


0.1S7E-03 


0. 17.jE-0.’i 


0. 302 E-03 


0.1 OS l>03 


0.7S5E-01 


O.olSE-Ol 


n.2.72E-03 


0.750E-0-1 


•10 


0.18SE-03 


o.i79E-o;f 


0.3 lOE-03 


0.20 3 E-03 


0 S73E-01 


0.509E-01 


0.303 E-03 


0.739E-01 


:>o 


D.lSlE-03 


U.IG7E-0;! 


O..331E-03 


0.200i:-03 


0.30SE-03 


0.4S7E-0 1 


0.411E-03 


0.7ME-01 


70 


0.20s E-03 


0.1()9E-0.3 


O.GO-lE-O;! 


0 217i;-()3 


0.106E-03 


0.479E-01 


0.62 IE- 03 


O.G3SE-(M 


90 


l).192E-03 


0 l l.')E-U.; 


0.7G')E-03 


0.222E-03 


0.772E-03 


0.442E-01 


0.178E-02 


0.688E-01 


110 


0.19SE-03 


0.1. 73 E-03 


0.3 so E-03 


0.220 E-03 


0.S74E-03 


0.333E-04 


0.18-1E-02 


0.533E-0-1 


130 


0.201 E-03 


0.155E-03 


0.3-11 E-03 


0.2 13 E-03 


0.S94E-03 


0 286E-04 


0.139E-02 


0.401E-01 


130 


0.201 E-03 


0 lGOE-0.3 


1 0.379E-03 


0 19.7E-03 


0.122E-02 


0.2.56E-04 


0.1.73E-02 


0 2S5E-04 



Table 3; 300 nodes with cquisi)accd arg'imcuts in (O.t); uniform weights. 





nodes in 


order of increasinc; argument. j 


nodes in random order. 




relative error in a,,, | 


error 


111 g,„ 


relative error m a,n 


error 


111 k-n 


It 


DD 


ED i 


DD 


UD 


DD 


UD 


DD 


UD 


0 


0.39SE-03 


0.39SE-03 


0.7 12E-U3 


0.7 12E-03 


0.1 2.7 E-03 


0.1 25 E-03 


O.llGE-03 


0.116E-03 


10 


0.457E-03 


0.398E-03 


0.18SE-02 


0.7.13 E-03 


0.125E-02 


O.lOlE-03 


0.106E-02 


0.121E-03 


20 


0.437E-03 


0 389 E-03 


0.187E-02 


0.720 E-03 


0.M4E-02 


0.112E-03 


0.146E-02 


0.122E-03 


30 


0.452E-03 


0.390E-03 


0.22SE-02 


0.7ME-03 


0.170E-02 


O.llOE-03 


0.163E-02 


0.131 E-03 


40 


0.4 79 E-03 


0.391 E-03 


0.265E-02 


0.694 E-03 


0.195E-02 


0.947E-0.1 


0.1 67 E-02 


0.113E-03 


.60 


0.508 E-03 


0 391E-0.3 


0.271 E-02 


0.6S0E-03 


0.215E-02 


0.886E-01 


0.1 65 E-02 


0.1 09 E-03 


70 


0.576E-03 


0.131E-0.3 


0 278E-02 


0 6-15E-03 


0.271 E-02 


0.796E-0I 


0.210E-02 


0.109E-0:! 


90 


0.533E-03 


0 4.70 E-03 


O.IGsE-02 


0.605E-03 


0.367E-02 


0.849E-0.L 


0.302E-02 


0.106E-03 


1 10 


0 558 E-03 


()..|37i:-0;t 


1) l!).3E-02 


0 752i:-0.3 


0.458E-02 


0.795E-04 


0.367 E-02 


0.102E-03 


130 


0.650E-03 


0 I2.:E-03 


' 0.2l0i:-02 


0 7()7E-():l 


0.72SE-02 


0.792E-0 1 


n.739E-02 


0.8 lOE-0 1 


1.70 


0.1 lOE-02 


0 .I9.1E-0.3 


j 0 376 E-02 


0 IG(il-:-03 


0.671 E-02 


0.548E-01 


0.728 E-02 


0.895E-01 
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Tabic -1: 300 nodes wilh landoin ar^maoiils in [0,2"); unifonn weights. 



i 


r(9ative oiTor m a,,, , 


error 


hi gm 


k 


1)D 


ID 


DD 


UD 


u 


u.i iii:-03 


1). 1 1 lL-03 


0.02.3 E-03 


U.025E-03 


\) 


0.50.si:-03 


U,ll)!Ji;-03 1 


0.0MTi:-03 


0.838E-03 


10 


0 J ISE-Ol 


l),lCUD-03 1 


; 0.31SE-01 


0.S38E-03 


•JO 


O.-JUC-Ol 


U.r."JD-03 


: 0.373E-01 


0.789E-03 


30 


0.-2*20E"01 


0.1()-lI--03 


0.32SE-01 


0.735E-03 


•10 


O.l ME-Ol 


U.117D-03 


0.131E-01 


0.826E-03 


50 


o.iS()i-:-oi 


0.1-.ME-0.3 


0.240E-01 


0. 6.54 E- 03 


70 


O.IOGE-Ol 


0.137E-03 


0.264 E-01 


0.G58E-03 


90 


O.JIOE-Ol 


0.rJ2E-03 


0.275E-01 


0.483E-03 


110 


0.317E-01 


0.137E-03 


0.820E-01 


O.lOOE-02 


130 


0.50-2E-01 


0 120E-03 


0.S32E-01 


0.508E-03 


150 


0.513E-01 


0.7^t!E-0l 


0.S6.5E-01 


0.360 E- 03 
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Table o: Dowiulaiing the FFT 



A' = 


1021 


rchiLi\o error iii a,„ 


error 


ill gm 


CPU seconds 


m 


Ic 


DD 


I'D 


DD 


UD 


DD 


UD 


1013 


11 


U.'l()lF-05 


0.213E-01 


0.2 It) E-05 


0.204E-01 


0.303L4-0 1 


0.148E-I-03 


093 


31 


0.502 E-05 


0.215E-01 


0.623E-05 


0.196E-01 


0.814E-I-01 


0.142E-1-03 


973 


51 


0.003E-05 


0.217E-01 


0.M2E-01 


0.193E-01 


0.132E-I-02 


0.136E4-03 


923 


101 


O.OOOE-05 


0.193E-01 


0.257E-01 


0.244E-01 


0.253E-i-02 


0.122E+03 


013 


111 


0.9-HE-05 


0 lOlE-01 


0.-2C5E-04 


0.242E-01 


0.276E-I-02 


0.119E4-03 


S93 


131 


O.lllE-O-l 


0.178E-01 


0 312E-01 


0.219E-01 


0.323E-r02 


0.114E+03 


.S73 


151 


0.120E-0-1 


0.17-1E-01 


0.373E-0-1 


0.220E-01 


0.36SE-I-02 


0.109E-(-03 


853 


171 


0.M3E-01 


0.172E-01 


0.-162E-04 


0.221E-01 


0.412E-h02 


0.104E-h03 


833 


101 


i 0.172E-01 


0.17-1E-01 


0.566E-04 


0.220E-01 


0.456E-^02 


0.990E-1-02 


>sl3 


211 


0.223E-0I 


0 171E-01 


0.721E-01 


0.241E-01 


0.497E+02 


0.941E-I-02 


713 


:U1 


0 177i:-0l 


II 1.33E-01 


0.131E-03 


0.169E-01 


0.692E-h02 


0.719E-H02 


013 


•111 


1 0.1G7E-03 


(I.103E-01 


0. 158E-03 


0.275E-01 


0.862E-802 


0.527E-H02 


513 


•511 


, l).2G7E-0.3 


l).709E-02 


0.527E-03 


0.174E-01 


O.lOlE-f-03 


0.365E4-02 


■113 


Oil 


0.28GE-03 


0.117E-02 


0 187E-03 


0.852E-02 


0.112E-f03 


0.233E+02 


313 


711 


0.3.37E-03 


0.322E-02 


0.551E-03 


0.7S1E-02 


0.122E+03 


0.132E+02 



A' = 
ni 


2018 

Ic 


relative error in i'i,„ j 

DD UD ' 


error 

DD 


in gm 
UD 


CPU : 
DD 


seconds 

UD 


2037 


11 


0.508E-05 


0.105E-00 


0.259E-05 


0.184E4-00 


0.391E-1-01 


0.390E-I-03 


1997 


51 


0.7-10E-05 


O.llOE-l-OO 


0.151 E-0 1 


0.235E-f-00 


0.172E4-02 


0.374 E-i-03 


• 1017 


101 


O.IOOE-IH 


i).102E-^00 


0 293E-04 


0.164E-i-00 


0.33GE-F02 


0 355E-f03 


; 1S07 


151 


0.125E-01 


O.O-tOE-Ol 


'I.12GE-01 


0.126E4-00 


0-49oE-l-0^2 


0.337E-f-03 


■ 1S'17 


201 


0.l5'JE-0 1 


0.9I8E-01 


0 (i 12E-01 


0.13lE-f00 


0.649E-f02 


0 319E-f03 


I 1797 


251 


('•2.37E-0-I 


O.ttlE-Ol 


1) t'l2E-01 


0,136E-i-00 


0.799E-h02 


0 301E-H03 


! 17*17 


:;oi 


0 2M1E-0 I 


Ht'(i5E-()l 


[ -.) 1 lCE-0.3 


0157E-800 


0.9K;E-+-02 


0 2^5E-i-0.3 


I 1007 


351 


o.:;5m:-o I 


0 t.3 lE-0l 


(.' 115E-03 


0.137E-I-00 


0.109E-f03 


0.2GSE-1-03 


1617 


•101 1 


0.5NlM)l 


0 sOOI'-Ol 


It 170E-03 


0.125E4-00 


0.123E-f-03 


0.252E+03 


1597 


•151 


0 749E-0 I 


0.7-ljE-Ol 


0 2 12E-03 


0133E-I-00 


0.13GE-I-03 


0.237E-h03 


1517 


501 


0.779E-0I 


0 729 E-0 1 


0.213E-0.3 


0.132E-1-00 


0.149E4-03 


0.2-22E-h03 


11*17 


001 


0 953 E-0 I 


0.G95E-01 


i 0.313E-0.3 


0.149E-I-00 


0.1 < ‘1 E-i-03 


0.19lE-f03 


13‘17 


701 


0.242E-03 


0.G28E-01 


0.39GE-03 


0.114E-i-00 


0.197E4-03 


0.1G7E4-03 


1M7 


001 


0.517E-03 


0 143E-01 


0.109E-02 


0.116E-I-00 


0.23SE-h03 


0.121E-f-03 


10*17 


1001 


0.731E-03 


0.381E-01 


1 0.113E-02 


0.640E-01 


0.25GE4-03 


0.100E-h03 


947 


1101 


0. 11 IE- 02 


0,'285E-01 


i 0.183E-02 


0.523E-01 


0.273E-h03 


0.815E-^02 
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